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Traditionally the history of remote sensing began during the First World War when aerial photography became a 
valuable reconnaissance tool. However, moving back more than a thousand years, the real pioneers of remote 
observation were probably the Nasca, a pre-Hispanic civilization living in southern Peru, between 100BC and 
700 AD. They used ‘earth observation’ as a mean of cultural expression drawing the geoglyphs (known as Nasca 
Lines) only visible from above. These drawings were made on flat desert surface of the Pampa by removing or 
clearing sand or stones, to create paths for ritual functions to please the gods and create harmonious relation¬ 
ships between man and environment. In this paper, the Nasca geo glyphs in Pampa de Atarco, are object of 
remote sensing based investigations with the twofold aim to identify and characterize them as well as to analyse 
and monitor their fragile state of conservation, threatened mainly by vandalism and off road vehicles. The 
approach herein proposed includes the integration and reuse of diverse remote sensing dataset, from multi- 
spectral satellite to Unmanned Aerial Vehicle (UAV) based LSAR data and close range photogrammetry. In 
particular, a multidate (2002-2013) very high resolution (VHR) optical satellite dataset has been processed in 
the spatial and temporal domain using textural indicators, including Skewness, Principal Component Analysis 
(PCA), and automatic classification tools which allowed us to enhance the visibility of disturbance features and 
to automatically extract them. The best results in terms of enhancement and automatic extraction capability of 
disturbance features have been obtained by Skewness. Moreover, the reuse of UAV L SAR-based correlation map, 
available free of charge from NASA, provided useful information on the state of disturbance from 2013 to 2015, 
widening the observation time window of the VHR satellite data set from 2002 to 2013. Finally, the integrated 
use of satellite VHR data with UAV-based photographs and DTMs, processed using structure from motion (SfM), 
allowed us to characterize, identify and reconstruct the relative chronological sequence of geoglyphs thus 
providing new insights and opening new perspectives for archaeological studies. 


1. Introduction 

Earth Observation (EO) technologies have multiple potentiality and 
numerous applications in the cultural heritage (CH) domain ranging 
from the knowledge improvement to the management (Lasaponara and 
Masini, 2008, 2017). In the last two decades the use of EO in CH has 
been strongly increasing thanks to the technological improvements of 
the sensors ever more useful for the study of the human past and an¬ 
cient landscapes for a large variety of environs including desert, tro¬ 
pical and Mediterranean ecosystems. In particular, the latest generation 
of satellite sensors opened up new frontiers and possibilities offering 
ever-closer and more comprehensive look at the earth's surface; pro¬ 
viding detailed imaging even for complex and fragile archaeological 
heritage as geoglyphs (Hesse, 2015; Masini et al., 2016b). 


As a whole, respect to the past, the archaeologists and CH decision 
makers are more aware of the benefits of EO technologies in terms of 
reduction of costs and time of archaeological investigations, and op¬ 
portunities to support strategies addressed to conservation and pre¬ 
servation of cultural assets thanks to the: 

(i) improvement in spectral and spatial resolution of sensors, evermore 
useful for site discovery and preventive archaeology; 

(ii) availability of multiscale, multitemporal and multi sensor data 
useful to investigate changes due to human activity in areas of 
cultural interest and monitoring risk in archaeological sites. 

Moreover, today these technologies are available at different costs, 
for different purposes and needs, and even with small budgets it is 
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possible to implement very effective solutions, but, they pose several 
issues related to the processing, integration, analysis and interpretation 
of data (Opitz and Herrmann, 2018). These issues are the challenges to 
face for transforming data into useful and reliable information and, 
therefore, they must be tackled by the scientific community in order to 
ensure an effective and reliable applicability of EO, as also expected in 
the Copernicus Programme initiatives (Casu et al., 2017). These require 
efforts also addressed to bridge EO technologies and archaeological 
communities to overcome the current limitations. In this context, sig¬ 
nificant contributions can be also achieved via case study demonstra¬ 
tion useful to support the development, application, assessment and 
improvement of reliable automatic data processing approaches and 
integrated methodologies of data analysis and interpretation. (Luo 
et al., 2019). 

In this paper, we focus on the analyses of Multitemporal/ 
Multisensor data, collected from satellite, UAV and ground, to obtain 
meaningful and detailed information about the changes that over time 
have affected the geoglyphs made by Nasca civilization in Southern 
Peru. They are only visible from above and are considered one of the 
most fragile and vulnerable world heritage sites (https://whc.u- 
nesco.org/en/list/700). This is because they were made: (i) removing 
and pulling the gravel inwards to the drawn lines, to obtain the figures 
in slight relief or (ii) through the removal of the gravel to obtain solid 
figures by the resulting contrast with its surroundings (Eitel, 2005); 
and, moreover, (iii) they sparsely spread up over huge areas with easy 
access, which make their destruction almost effortless and their pro¬ 
tection very challenging (Masini et al., 2016b). These characteristics 
highlight the challenges and urgent needs related to their protection 
and preservation issues being that Nasca geoglyphs provide a unique 
and invaluable source of information about Nasca civilization that 
flourished in the distant past. For these reasons, specific surveillance 
and monitoring strategies must be adopted to safeguard them. 

The methodology herein devised is based on the use of hetero¬ 
geneous information from multiple EO data sources ad hoc processed to 
improve knowledge and set up monitoring methodologies of geoglyphs 
and their surrounding landscape. Knowledge improvement is a man¬ 
datory step for the preservation especially for the specific character¬ 
istics of geoglyphs and Nasca Lines considering the difficulty in re¬ 
cognizing and detecting them in particular for those located in Pampa 
de Atarco. In fact, the specific drawing techniques (see section 2.3) of 
these geoglyphs make them very difficult to see not only from on 
grotmd, but also from above (from airplane, drone and satellite). For 
this reason, Pampa de Atarco (see section 2) has been herein selected as 
test site, to develop a suitable methodology, based on remote sensing 
(the only tool applicable for geoglyphs!), to enhance the visibility of 
geoglyphs and to capture and extract changes that adversely affecting 
them. 

Our effort is a contribution to the definition of the best practices for 
the production of critical information necessary for improving the 
knowledge and set up a reliable and systematic monitoring. This is 
pursued by: (i) exploiting historical and updated data available from 
remote sensing technologies from satellite, manned and unmanned 
aerial vehicles (UAVs), (ii) defining ad hoc data processing to enhance, 
detect and characterize geoglyphs features, (iii) using key statistical 
indicators to capture and quantitatively evaluate the impact of changes, 
(iv) the validation of the outputs from satellite data by using close range 
data. In particular, the use of statistical analysis (PCA and Skewness) 
along with unsupervised classification allowed us not only to identify 
and delimiting disturbance, but also to extract in automatic way sui¬ 
table information useful to face the ongoing threats and facilitate the 
operational activities and preservation strategies. 

The methodology herein proposed can be promptly re-applied to 
other cultural heritage sites with particular reference to those char¬ 
acterized by geoglyphs which were a cultural phenomenon common to 
several civilizations in different historical periods and geographical 
regions. In particular, the drawing of figures only visible from above, 


was very common in Americas, see, for example, the serpent-mound in 
Ohio, the ornithomorphic figures in Wisconsin (Orefici, 2009), the 
geoglyphs in the Atacama desert of Chile (Briones, 2006), channels and 
embankments depicting geometric figures in the Acre state (Brazil) 
(Saunaluoma and Schaan, 2012), and, finally, the ring ditches in the 
Bolivian Amazon (Erickson, 2010). In Peru, the geoglyphs were drawn 
usually on the desert coast, where the sandy loam soil, covered by al¬ 
luvium, was more suitable for their conservation, see, for example, the 
Condor de Oyotun and the anthropomorphic figure in Pampa de Cana 
Cruz in the Lambayeque region (Alva and Meneses Alva, 1984), the 
geoglyphs of Canto Grande, near Lima (Rosello et al., 1985), and, of 
course, the well known Lines and Geoglyphs of Nasca and Palpa, listed 
in World Heritage since 1994 (https://whc.unesco.org/en/list/700). 

2. Study area and data set 

2.1. The geoglyphs of Nasca: technical aspects and cultural meaning 

The Nasca Lines are one of the most impressive cultural heritages in 
the world, dated from 200BC to 600 AD and created by Paracas and 
Nasca civilizations considered among the richest, most refined, and 
original cultures of the pre-Hispanic Andean world (Silverman, 1990; 
Masini et al., 2016). These geoglyphs were obtained using three dif¬ 
ferent engraving techniques by: (i) removal of gravels and sand to 
evidence the underlay lighter dust, (ii) adding dark colour pebbles and 
sand grains along the lines, (iii) creating slight micro-relief by scraping 
sand and gravels and adding some darker gravel along the edges of the 
figures (Eitel, 2005; Masini et al., 2016a). 

The impressive geoglyphs shapes, huge sizes (lines longer more than 
1 km and biomorphic figures reaching up to 300 m; see Fig. lbcd), and, 
functions have attracted an extraordinary scientific interest. In the last 
century, diverse assumptions and numerous theories have been postu¬ 
lated among them the main is the astronomical one, which considers 
the Nasca lines as a “great book” of astronomy with the function of a 
calendar (Kosok and Reiche, 1949; Aveni, 1986). The latest theories are 
based on the spatial/cultural/religious relationship between (i) the 
mountains and sacred places (Reinhardt, 1988); or (ii) between geo¬ 
glyphs and settlements (Reindel and Wagner, 2009), and Cahuachi 
ceremonial center (Silvermann, 1990; Orefici, 2009; Masini et al., 
2016a). According to these latest theories, the geoglyphs were sacred 
spaces used as “venues” for civic ceremonial activities and pilgrimage 
paths. 

As a whole, the phenomenon of the geoglyphs has its origin in a 
profound religious and cultural substratum, strongly influenced by the 
extreme conditions of the climate and environment where the Nasca 
used to live. Therefore, the geoglyphs were a cultural expression of the 
need of adaption, a request of help to the gods, an attempt to create 
harmonious relationships between man and environment, earth and sky 
made by ceremonial activities performed in the geoglyphs sacred space 
(Masini et al., 2016a,Masini et al., 2016). Geoglyphs were made to 
please and to be seen by Gods, therefore, only visible from above. 

2.2. Protection issues 

Due to their specific characteristics, geoglyphs are very fragile be¬ 
cause they can be severely damaged or completely destroyed simply 
stepping on them. Therefore, the protection of this vast area is a 
mandatory issue to reduce threats mainly due to tourism and vandalism 
(UNESCO, 1994). For these reasons, Peruvian Law No. 24047/19, put 
the National Institute for Culture in charge to protect these geoglyphs, 
the access was forbidden and the planning for urban and rural devel¬ 
opment limited. But, despite this, unfortunately, over the years, nu¬ 
merous disturbances have left indelible traces on the landscape as, for 
example, tracks made by vehicles or feet, including those linked to the 
Dakar auto rallies of 2012 and 2013, and other damage (Kozak, 2014). 

In 1993, the Lines and geoglyphs of the Pampas of Nasca spreading 
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Fig. 1. . (a) Landsat TM acquired in 2008 of Rio Grande drainage basin, in Southern Peru, including the Pampa de Jumana characterized by the most famous 
geoglyphs (among which those showed in Fig. lb and c), Pampa de Atarco which is the area under investigation of this paper, the Ceremonial center of Cahuachi, and 
the tributaries of the river, among which Rio Nasca. (b-c) Aerial images (© Nicola Masini) of the famous hummingbird and a spyraliform motif located in Pampa de 
Jumana. (d) The hummingbird seen from the ground (courtesy by Giuseppe Orefici). 


over more than 700 square km2 along low foothills and desert terrain 
were declared Archaeological Reserve Area and later on in 1994 in¬ 
scribed on the World Heritage List. According to the 2013 UNESCO 
report, the factors affecting the geoglyphs are: i) illegal farming and 
mining activities; ii) vehicular tracks; iii) flooding; iv) lack of systematic 


monitoring; v) insufficient air-traffic security measures; vi) lack of a 
management plan. 

On the January 16, 2015, a Resolution of Peruvian Ministry of 
Culture (No. 019-2015-MC) approved the “Management System for 
Cultural Heritage Nasca and Palpa Territory” plan that covered not only 
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Fig. 2. (a), (b) Satellite based DEM and orthophoto of Pampa de Atarco geoglyphs, Nasca River and Cahuachi with identification of the area under investigation, 
indicated with green box. Red and yellow polygons indicate the areas of aerial and drone surveys, respectively, (c) Orthophoto obtained processing aerial images by 
SfM. (d) Detail of aerial orthophoto whose location is indicated with blue box in Fig. 2d. (e) Detail of UAV-based DEM whose location is indicated with red box in 
Fig. 3c. (f-g) Profiles X-X is related to a trapezoidal geoglyph, whereas profile Y-Y is related to a disturbance track (see also Fig. 2d); a and b, in Fig. 2f, refer to the 
edges of the trapezoidal geoglyph, c refers to the track (see also Fig. 2d and e). In Fig. 2d green box indicates the area that will be showed in Fig. 10. 
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the UNESCO mapped area but also a larger territory (over 5627 km2) 
which also included the geoglyphs of Pampa de Atarco object of ana¬ 
lysis in this paper. 

2.3. Study area: pampa de atarco geoglyphs 

Pampa de Atarco, located at 2-4 km south of the pyramids of the 
ceremonial center of Cahuachi, at 5-7 km south of the Nasca river (see 
Fig. la, 2a-c), fits into the hollow of the Rio Grande de Nasca river 
systems, tributaries of the Rio Grande, dating from the Tertiary and 
Quaternary periods. It lies on an important impermeable base made up 
of metamorphic, sedimentary, and intrusive rocks (Delle Rose, 2016; 
Orefici and Lancho Rojas, 2016; Schreiber& Lancho Rojas, 2009). 

Previous studies (Masini et al., 2016a) highlighted a spatial and 
functional relationship between the geoglyphs of Atarco and the tem¬ 
ples of Cahuachi. Both Geoglyphs and temples were sacred places used 
by the Nasca for religious ceremonies and other liturgical forms of re¬ 
ligion; hence the cultural importance of the geoglyphs of Atarco. 

As regards the morphological characteristics, the presence of geo¬ 
metric figures is prevailing as straight lines, trapezoidal, rectangular, 
spyraliform, and meandering motifs. From the execution point of view, 
all the three different techniques, already described in 2.1, are present. 
In Pampa de Atarco, the use of sand grain, pebbles and gravel of small 
dimension (gravel size ranging from 8 to 20 mm, pebble size from 70 to 
120 mm; see Fig. 3b) make the geoglyphs very difficult to see on ground 
(Fig. 3a) and, in not optimal illumination conditions even from air craft 
(Fig. 3c) or drone. This makes necessary the use of diverse remote 
sensing dataset and processing techniques aimed at enhancing edges 
and texture of the geoglyphs. 

2.4. State of the art of remote sensing of Nasca geoglyphs 

Despite the great interest in studying geoglyphs using observation 
tools from above, up to now, a few attempts have been only made to 
evaluate the performance of EO for the monitoring, protection, and 
preservation of geoglyphs of Nasca, and Palpa. In the last decade, a few 
but significant studies have been conducted using SAR earth observa¬ 
tion technologies. 

Ruescas et al. (2010) used the radar coherence to assess changes of 
the Nasca Lines from 1997 to 2004. In particular, they used two pairs of 
satellite SAR data: the first acquired by ERS-2 SAR from 1997 and 1999; 
the second one by ENVISAT-ASAR from 2003 and 2004. Several dif¬ 
ferent decorrelation factors contributing to a loss of coherency in a 
radar pair can be distinguished, and these include the temporal change 
in the ground properties and nature between the two satellite passes 
(Ruescas et al., 2010). 

Chapman et al., (2015) evaluated the potential of NASA/JPL UAV 
SAR for examining the condition of the Nasca geoglyphs. In particular, 
L-band data were collected almost exactly two years apart, on 19 March 
2013 and on 21 March 2015. In particular, they used the repeat-pass 
interferometric pair to detect disturbance affecting the Nasca lines, 
including erosion phenomena. Following Chapman et al., 2015, the 
same author group (Comer et al., 2017) integrated the L-band from 
UAV with Sentinel-1 C-band to detect landscape disturbance threa¬ 
tening Nasca and Palpa Lines. Cigna and Tapete (2018) used COSMO- 
SkyMed InSAR for tracking some human-induced disturbance close the 
famous hummingbird. 

As regards optical data, Lambers (2006) used aerial photo- 
grammetry for the documentation and morphological characterization 
of the Palpa geoglyphs. Hesse (2015) investigated a Landsat TM time 
series (2009-2013) to identify disturbances affecting the Lines. He 
pointed out that the most extensive surface disturbances occurred be¬ 
tween 06 January 2012 and 22 January 2012, due to the intense off¬ 
road vehicle linked to 2012 Dakar rally which passed close to the area. 
Finally, very high resolution optical satellite imagery have been used by 
Masini et al. (2016a,b) for mapping the Ceremonial centre of Cahuachi 


and the geoglyphs of Pampa de Atarco, with the twofold aim to provide 
new information on the function of these geoglyphs and to map dis¬ 
turbance caused by off-road vehicles and urban sprawl. The approach 
was simply based on the visual reconnaissance of both geoglyphs and 
disturbances affecting them. 

In this paper, we focused on the challenging issue of setting up a 
multiscale and multisensor methodological approach (i) to facilitate the 
visual identification of geoglyphs of Atarco, which are much more 
subtle and difficult to be recognized (see Fig. 3) compared to those of 
Palpa and Nasca, and (ii) to detect changes linked to anthropic dis¬ 
turbances, evaluating the opportunities to extract them in an automatic 
way. 

3. Methodology 

The approach herein proposed is based on the identification of key 
statistical indicators useful (i) to enhance and characterize geoglyphs, 
and (ii) to automatically detect and extract changes linked to dis¬ 
turbance. To this aim, we adopted an integrated approach, based on 
diverse remote sensing technologies and analyses, including close range 
RS data. On the basis of the availability of data, we used a multi¬ 
temporal (2002, 2005, 2011 and 2013) optical VHR data set and UAV 
L-band SAR survey acquired in 2013 and 2015 and processed by NASA. 
Finally, aerial and drone surveys were conducted in September 2015 in 
a significant subset of the investigated area. The analyses were con¬ 
ducted using the optical satellite data listed in Table 1. 

After the preprocessing (calibration, atmospheric correction, or¬ 
thorectification), the satellite optical VHR images were processed using 
multi-temporal Principal Component Analysis (PCA) and texture to 
enhance, speed, and improve the identification of both geoglyphs and 
disturbances. In particular, the information on the texture is herein 
performed in both spatial and temporal domains. Mathematical pro¬ 
cedures used to extract texture features, can be categorized into three 
general approaches: (i) statistical analysis (as for example PCA), (ii) 
model-based methods, and (iii) signal processing methods, including 
Fourier transforms, convolution filters, co-occurrence matrix, spatial 
autocorrelation, fractals, etc. In this paper, the investigations were 
performed as follow: 

(a) Four VHR satellite optical images were processed using multi¬ 
temporal analysis based on: (i) PCA (ii) textural indicators, in¬ 
cluding Skewness, and (iii) unsupervised classification of PCs and 
Skewness; 

(b) Re-use of two date acquisitions of L-SAR UAV as processed by NASA 
to obtain the interferometric correlation (Chapman et al., 2015). 

(c) Processing of RGB images taken from air flight and drone using 
Structure from Motion (SfM), to obtain orthophotos, digital eleva¬ 
tion models (DEMs), and further enhancement of DEMS-micro-to- 
pography using Simple Local Relief Model (SLRM) and Sky View 
Factor (SVF); 

The processing of VHR satellite data was only based on the pan¬ 
chromatic scenes because, according to previous investigations of the 
same author group (Masini et al., 2016a), in this desert area, the mul- 
tispectral channels do not provide significant additional information 
compared to panchromatic and, therefore, we preferred to exploit the 
higher spatial resolution. 

3.1. Processing the optical very high resolution (VHR) multi-date satellite 
imagery 

3.1.1. Principal components analysis and skewness 

The need to identify and characterize changes mainly due to human 
induced disturbance (from footpaths to vehicle tracks) makes necessary 
to perform multitemporal analyses. To this aim, PCA and textural 
analyses, based on Skewness, have been used. 
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Fig. 3. All ground operations (2a), including drone driving (2e), have been carried out with all possible precautions such as the use of particular shoes, GPS with 
satellite images and maps in order to avoid to damage geoglyphs stepping on them. These images show how fragile are the geoglyphs and how difficult to identify 
them on ground (2b, 2c, 2d) and by drone (3f). The ground images (2a-2e) are by Rosa Lasaponara; The oblique aerial image taken from UAV is by Antonio Pecci. 


The PCA (see, for example, Richards and Jia, 2006; Neeti and 
Eastman, 2014) is a mathematical transformation that is useful and 
used to remove redundancy in a given multi-band (multispectral or 
multitemporal) data set through the creation of a new image data set 
with fewer uncorrelated components. This is obtained by translating 
and/or rotating the axes of the original feature space, using a linear 


transformation to decorrelate multivariate data and represent them in a 
new component space without correlation. 

To do this, the process computes the covariance (or variance) matrix 
(S) among all input time series, the eigenvalues, eigenvectors and 
loadings using the formulas shown in the Appendix. 

The most common application of PCA is based on the use of 
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Table 1 

Satellite optical VHR data set: acquisition time, satellite sensors, ground sample distance (GSD). 


Acquisition date 

Satellite sensor 

Panchromatic GSD (m) 

Multispectral GSD (m) 

Time 

Other 

16.09.2002 

QuickBird-2 

0.619 

2.48 

15:17 


25.03.2005 

QuickBird-2 

0.634 

2.54 

15:29 


28.02.2011 

Stereo GeoEye-1 * 

0.50 

2.00 

15:30 

DEM 

13.03.2013 

Pleiades IB 

0.50 

2.00 

15:20 



covariance matrix, even if both variance and covariance measure the 
“spread” of a dataset around the mean. A series of new image layers 
(eigenvectors generally called eigenchannels or components) are ob¬ 
tained from the eigenvalues of the covariance matrix. 

PCA allows us to highlight the area (pixels in this case) affected by 
changes, due to the low correlation associated (in multi-temporal da¬ 
taset) to the pixels that change substantially. The major portion of the 
variance in a multi-temporal image data set is associated with non 
variable surface characteristics and will be mainly in the first compo¬ 
nent; whereas, pixels/areas of localized change will be enhanced in 
later components. To correctly interpret the results obtained from the 
PCA additional information (as loadings, ground survey, etc.) are re¬ 
quired being that the eigenvectors are extracted from the series itself 
and, therefore, they cannot have a general and universal meaning. 

Along with PCA, we also considered the Skewness (Balakrishnan 
and Scarpa (2012), applied in the temporal domain, in order to detect 
changes. Skewness measures the degree of distortion from the sym¬ 
metrical normal distribution informing on how much and in which 
direction, the values deviate from the mean, median and mode. The 
Skewness can be determined by how these quantities are related to each 
other, and for this purpose there are several formulations that are also 
available in image processing tools (https://www.itl.nist.gov/div898/ 
handbook/eda/section3/eda35b.htm). The most common skewness 
measurement is Fisher-Pearson (or Pearson median skewness) ex¬ 
pressed by the formulas reported in the Appendix. 

3.1.2. Unsupervised classification and segmentation 

To automatically extract changes, outputs from PCA and Skewness 
were further processed using unsupervised classification selected be¬ 
cause: (i) it is an automatic process, (ii) classes do not have to be de¬ 
fined a priori; and, moreover, (iii) unknown features/classes may be 
identified and categorized. 

Unsupervised classifications can be performed using different ap¬ 
proaches/algorithms, as: (i) K-means clustering, and (ii) ISODATA 
(Iterative Self-Organizing Data Analysis Technique) (Ball and Hall, 
1965). In this study, ISODATA was applied because more flexible 
compared to k-means which a priori assumes the number of clusters. 

The classification is followed by a segmentation step to reduce the 
complexity of the classification outputs and facilitate the image inter¬ 
pretation. The segmentation process requires the following three 
parameters: (i) Colour, (ii) Scale, and (iii) Form (Shape). 

The first parameter balances the homogeneity of segment colour 
and shape: value “one” will result in very fractal segments (with a low 
standard deviation), whereas a zero value in very compact segments 
(with higher colour heterogeneity). 

The scale parameter is an abstract value without any direct corre¬ 
lation with the object size, and it depends, rather, on the heterogeneity 
of the data, performed setting: i) the minimum number of pixels for 
building a segment (in our case set equal to 50; ii) the number neigh¬ 
bouring pixels which determines the separability/connectivity of the 
segments (set equal to 4). 

Finally, the form parameter controls the form features balancing the 
criteria for border smoothness and object compactness. 

3.1.3. Spatial texture analysis 

Texture analyses have been also performed in the spatial domain to 


enhance edges of geoglyphs features and, therefore, to facilitate their 
identification from above. Texture is not an absolute parameter, but, it 
is function of spatial and radiometric resolution and according to the 
given scale can be fine, rough, rippled, smooth, and irregular (Irons and 
Petersen, 1981). In this paper, texture analysis has been performed 
using co-occurrence measures and the relative outputs are based on 
several indicators such as mean, variance, dissimilarity, entropy, etc 
(see equations (7) to 12) shown in the Appendix) to quantitatively 
evaluate the textural parameters. The texture categories herein used are 
based on distinguishing features obtained from the Mean and Standard 
Deviation which are the local mean and standard deviation values ob¬ 
tained for the selected size of the processing window, as well as from 
the Variance, Correlation, Contrast, Dissimilarity, Homogeneity, and 
Entropy computed as shown in the Appendix. These indicators are all 
computed on a single scene and then compared with the output ob¬ 
tained over time from all the available multi-temporal data sets (ac¬ 
quired on 2002, 2005, 2011, 2013; see Table 1). Outputs form texture 
analysis have been classified using the approach describes in section 
3.1.2. 


3.2. Reuse of UAV SAR interferometric coherence 

The UAVSAR data were processed by NASA (Chapman et al., 2015; 
Comer et al., 2017) and available online free of charge (https://uav- 
sar.jpl.nasa.gov/cgi-bin/data.pl). Among the products, we used the 
correlation image computed from the SLC data, to identify the changes 
occurred between the two years of acquisition 2013-2015 in order to 
integrate those mapped from the VHR satellite between 2002 and 2013. 
The correlation was computed as follows using formula: 

A _ + SNR- 1 )) 

i s i..i 2 -zr=ii s 2i 2 [i3] 

where. 

N is the number of pixels in the coherence window calculated, 

Si„ and S 2n indicate the complex value of master and slave image, 
respectively, 

indicates conjugate operator, and 
SNR is the signal-to-noise ratio 

The correlation y can vary between 0 (complete decorrelation or, 
equivalently, no correlation between the two images) and 1 (no dec¬ 
orrelation or, equivalently, complete correlation between the images). 
The need to normalize y by the SNR is due to the fact that some areas, 
shadowed by topographic features, may have a correlation greater than 
1. 

Some authors have reported promising results from investigations 
on the relationship between InSAR coherence (correlation) and changes 
in superficial cover and/or other characteristics (Chen et al., 2016). It is 
expected that the use of interferometric may enhance changes, thus 
improving the identification and characterization of archaeological 
remains. 
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3.3. Manned and unmanned aerial image acquisition and processing 

Close range acquisitions were performed in 2015 taking photo¬ 
graphs using a RGB camera from an ultra light aircraft (see Table 2). 
The aim was to identify subtle geoglyph features and characterize the 
superposition of drawings in order to understand the relative sequence 
by the interpretation of oblique photos, DEMs and orthophotos. 

The surveyed area was the center of the meandering figure. Some 
detailed images were acquired using a drone Dji Phantom Vision 2, 
mounting a 14 megapixels RGB camera. 

Both manned and unmanned aerial images have been processed 
using SfM algorithm in Agisoft Photoscan software ver. 1.4. 

To appreciate the microtopography with centrimetric detail for the 
interpretation of geoglyphs and disturbance features, the DEMs ob¬ 
tained from both manned and unmanned images were enhanced using a 
number of visualization techniques, including hill shading, simple local 
relief model (SLRM) and Sky View Factor (SVF). 

From the computational point of view, SLRM (Hesse, 2010) ap¬ 
proach is based on the: (i) smoothing of DEM made applying low pass 
filter, (ii) its subtraction to the initial DEM, (iii) calculation of the zero 
meter contours from the difference model to obtain break lines, as well 
as the intersection of the break lines with the DEM. SLRM computation 
is based on the use of a kernel within computing the statistical para¬ 
meters and, obviously, the results will strongly depend on the kernel 
selection. 

SVF (Zaksek et al., 2011) quantifies ‘the portion of the sky visible 
from a certain point’ within a certain radius, in other words large SVF 
denotes a large portion of the sky visible and the size of the observed 
area (defined by the chosen radius) impacts on the result. Small radius 
is required (e.g., 10-15 m) to highlight subtle and small-scale micro¬ 
relief as expected in the geoglyph area. In particular, SVF delineates 
mainly concave features and maintain the slope well visible also pre¬ 
serving the perception of the general topography. 

4. Results and discussion 

A visual inspection of the multitemporal satellite images clearly put 
in evidence that these data, without a specific data processing, do not 
provide the details necessary to (i) analyse the morphological features 
of the geoglyphs and (ii) identify changes affecting them (with parti¬ 
cular reference to anthropic disturbances). As an example, Fig. 4(a-d) 
shows a detail of a meandering geoglyph, thought to be linked to cer¬ 
emonial activities performed by the Nasca people, walking in proces¬ 
sion, singing and praying on the ground drawing (Masini et al., 2016a; 
Orefici, 2009). This scene, from now on, will be assumed as reference 
for the following figures. It is characterized by the presence of geo¬ 
glyphs shaped as lines (in particular a trapezoid) and curvilinear ele¬ 
ments, obtained by the removal and addition of stone material to create 
micro-relief (see Fig. 2 in section 2.3 and also details from UAV derived 
DEM in Fig. 3d-f). 

As regards the disturbance features, only the big tracks caused by 
off-road vehicles are visible, as indicated in Fig. 4 by coloured arrows 
and letters dl, d2, d3, d4, and d5. They have a wide that ranges from 4 
to 16 m. As a whole, we can argue that the visual comparison of the 
panchromatic scenes shows the biggest changes and a different fre- 
quentation of the tracks, in one case continuous (see d2 in Fig. 4), in 
other cases more sporadic. 

4.1. Multitemporal analysis for disturbance identification and automatic 
extraction 

Some statistical computations (using ENVI tools) have been per¬ 
formed to enhance geoglyphs, and changes affecting them, and to 
overcome the subjectivity of the visual comparative analyses. In parti¬ 
cular, three strategies have been adopted to extract information and 
translate data into useful, operational, and meaningful information. In 
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Fig. 4. (a) 2011 GeoEye panchromatic image. The green box indicates the whole area under investigation. The red box is the subset herein discussed and showed as 
multitemporal date set in b-e. In particular, (b) 2002, (c) 2005, (d) 2011, (e) 2013. The coloured arrows and the letters dl-d5 denote disturbance tracks caused by off¬ 
road vehicles. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) 


particular, we used: 

i) PCA and Skewness to identify and characterize disturbance features 
from 2002 to 2013 

ii) Unsupervised classification and segmentation to automatically 
identify and map changes 

iii) SAR correlation map for estimating the changes occurred from 2013 
to 2015 

Outputs from PCA are in Fig. 5 which shows the four principal 
components as obtained from the multitemporal panchromatic satellite 
dataset. 

The first component (PCI) accounts for the majority of the total 
dataset variance exhibiting strong correlation values among the four 


input scenes. 

As expected, being that PCI is an integral/average value, and con¬ 
sidered that geoglyphs are quite stable over time (for more than one 
thousand years ..!), it enhances geoglyphs and potential permanent 
disturbance, as in the case of the central black track (indicated with 
orange arrows and letter d2) caused by off-road vehicle, south-north 
oriented, and visible since 2002 (see Fig. 4). It has a wide that ranges 
from 3 to 5 m and is present in all the scenes (from 2002 to 2013, see 
Fig. 4a-d) thus indicating that it has been used continuously. On the 
basis of the previous consideration, we can point out that PCI well 
enhances the geoglyph edges and maintains the traces of disturbances 
permanently present in the whole dataset (as in the case of the feature 
denoted as d2 in Fig. 4a). 

Later components (see Fig. 5b-d) account for the temporal 
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Fig. 5. The four principal components as obtained from the PCA applied to the multitemporal panchromatic satellite dataset. The arrows indicate the main dis¬ 
turbance features. 


Table 3 

Summarizes the reconnaissance of all the disturbance features from the satellite 
scenes and PCA components. 


Disturbance features (vehicle tracks) 

Subtle disturbance features 
(footpaths or motorcycle 
tracks?) 

DATA dl 

d2 

d3 

d4 

d5 

x y z k w 

2002 

X 

X 


X 


2005 

X 

X 


X 


2011 x 

X 


X 



2013 x 

X 





PCI 

X 





PC2 x 

X 

X 

X 


X 

PC3 

X 

X 


X 

X X 

PC4 x 

X 


X 


X X 


variability of the image multi-temporal dataset and well emphasize 
changes even small ones, as in the case of the tracks clearly highlighted 
in PC2, PC3 and PC4. Compared to the panchromatic scenes, some 
tracks are much more visible, as those indicated by purple, cyano, red 
and green arrows and letters dl, d3-d4 in Fig. 4b-d. As a whole, the 
later PCA components well evidence all the disturbance features, in¬ 
cluding the ‘permanent’ one (d2), due to the fact that over the years 
vehicles do not follow exactly the same path, causing track changes in 
both shape and width well revealed by the statistical de-correlation. As 
regards the other non-permanent disturbance features (dl, d3-d5), PC2 
provides better results compared to PC3 and PC4 (see Table 3). 

Ever more important is the fact that the later components put in 
evidence disturbances not visible from the four panchromatic images 
(please compare Fig. 5b-d with 4a-d). In particular, they are subtle 
features, probably linked to footpaths or motorcycle tracks, indicated 
with black arrows and letters x, y, z, k, w (see also Table 3). Moreover, 
we notice the ‘complementarity’ of PC2, PC3 and PC4 in evidencing the 
above said disturbance traces. In particular, the feature ‘y 1 is only visible 
from PC2, ‘x 1 and ‘w’ are observed only in PC3; and finally, the other 
two indicated with ‘z 1 and ‘k 1 compare only in PC4 (see Table 3). 

As a whole, respect to the single panchromatic scenes (see Fig. 4) 
PCA provides a more complete and detailed survey of disturbance 
features. In particular, it evidences all the five disturbance features 


observed also from the panchromatic satellite dataset, plus five subtle 
disturbance tracks (indicated with black arrows and letters x, y, z, w, 
and k). 

The four PC components were classified and segmented to auto¬ 
matically extract the changes. The results are not fully satisfactory. In 
particular, for each component a few disturbance tracks along with 
erosion are extracted. Fig. 6a shows the automatic extraction of PC2 
which exhibits only three disturbance tracks and other changes linked 
to areas involved in erosion processes. 

For this reason, further analyses have been performed using multi¬ 
temporal texture analyses computed on the whole dataset such as the 
Standard deviation and the Skewness. The best results have been ob¬ 
tained from Skewness as showed in Fig. 6b-d and Table 4. In particular, 
we can notice that Skewness enables the visual identification of all the 
five ‘strong’ disturbance features caused by off-road vehicles (indicated 
with coloured arrows), and three of the five subtle one. On the eight 
anthropogenic disturbance recognized by visual inspection, seven are 
automatically extracted along with areas affected by changes probably 
due erosion processes. 

From the rationale point of view both PCA and skewness indicators 
measure how much and in which direction, the values deviate from the 
mean. In particular, Skewness produces one image output that informs 
us on how much each pixel value deviate from the mean, whereas, PCA 
produces as output various components (whose max number is the 
number of the input bands), which enable us to categorize how much 
each pixel value deviate from the mean. Each successive component 
contains less of the total dataset variance and using the loadings (ob¬ 
tained from formula 4) PCA also provides information related to when 
the change occurred (namely in each input band in our case multi¬ 
temporal band). From the rationale point of view, we propose the joint 
use of PCA and Skewness, in order to (i) immediately capture the 
change even subtle using the skewness (ii) maintaining from PCA the 
information related to how much each pixel value deviate from the 
mean and when this occurred. Moreover, the Skewness can enable us to 
automatically categorize all the changes in an easier way compared to 
the diverse PCA components, that can be characterized by clusters 
which present low internal heterogeneity and, in turn, can be more 
complex to automatically categorize changes using an unsupervised 
classification. 
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Fig. 6. (a) Automatic extraction of PC2 which exhibits three linear disturbance features, two caused by off-road vehicles and two, smaller and subtle, due to footpaths 
or motorcycle tracks, (b), (c) and (d) are related to Skewness, its unsupervised classification and segmentation, respectively. The Skewness allowed us to identify 
eight disturbance features, among the six are automatically extracted. 


As a whole, we can point out that: (i) the use of both PCA and 
Skewness enabled us to well emphasize changes over time with parti¬ 
cular reference to those linked to anthropogenic disturbances. PCA 
components provide more details, related to the number of the change 
features, their strength of de-correlation and the time when occurred, 
but, the automatic extraction is really poor being that weak categor¬ 
izations are generally provided by all the PCs (except PC2). On the 
other side. Skewness outputs enhance a fewer number of change fea¬ 
tures compared to PCA components (eight over ten compared to) but it 
is suitable to automatically categorize most of changes numerically 
seven over eight. 

For estimating the changes occurred after 2013 we re-used the 
correlation map computed by NASA using the L-band data acquired on 
19 March 2013 and on 21 March 2015 from UAV SAR platform. To 
facilitate the interpretation of the correlation map we overlaid it (see 
Fig. 7) on the map of geoglyphs obtained as explained in section 4.2. In 
Fig. 7, the green colour denotes areas with higher correlation, and 
therefore not affected by significant changes. Whereas, yellow and or¬ 
ange colour characterize those areas with lower correlation values 
ranging from 0.5 to 0.7, thus denoting significant variations. In parti¬ 
cular, we notice two linear patterns which follow two disturbance 
tracks (indicated by black arrows, and d2 and d5) already identified 
from the multitemporal dataset and relative outputs (PCA, textural 
analysis; see Figs. 4-6)). Such correlation values indicate that some 

Table 4 

Summarizes the rate of visibility and automatic extraction of the disturbance features from PC2 and Skewness. 


Disturbance features (vehicle tracks) Subtle disturbance features (footpaths or motorcycle tracks ?) 


DATA 

dl 

d2 

d3 

d4 

d5 

X 

y 

z 

k 

w 

PC2-based reconaissance 

X 

X 

X 

X 



X 




PC2 

X 



X 



X 




Autom. extract 











Skewness based 

X 

X 

X 

X 

X 

X 

X 



X 

Reconaissance 











Skewness 

X 

X 

X 

X 

X 

X 

X 




Autom. extract 












changes occurred between 2013 and 2015. Consequently the dis¬ 
turbance activity continued after the last Dakar rally. Moreover, scat¬ 
tered disturbance features due probably to footpaths (pd) characterize 
the trapezoid and some parts of the meandering geoglyph. At northeast 
and southwest of the scene are indicated with letter ‘e’ two areas af¬ 
fected by changes reasonably due to erosion processes. 

As a whole, based on the results of the processing of multitemporal 
panchromatic dataset from 2002 to 2013 (see Figs. 4-6) and the output 
from the interpretation of SAR correlation map (2013 and 2015), we 
can conclude that respect to the ten disturbance tracks, five reasonably 
caused by off-road vehicles (see dl-d5 in Figs. 4-6) and 5 linked to 
motorcycle tracks or footpaths, only two, vehicle tracks (d2 and d5) 
were still in use after 2013. Local disturbance features are probably due 
to wind erosion and vandalism. 

4.2. Enhancement and characterization of geoglyph features 

In respect to the geoglyph features, considering the contribution of 
PCA, only the first component enabled us to enhance them. In parti¬ 
cular, as expected, their visibility is improved by PCI compared to the 
single scenes (compare Fig. 8a with Fig. 4a-d). In fact, focusing on the 
single scenes we notice a diverse visibility of the edges of geoglyphs, 
from one scene to another. For example, the trapezoid is more visible in 
2002 and 2005; whereas, some of the curvilinear elements are more 
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Fig. 7. UAV-SAR based correlation with geoglyphs features mapped as detailed in section 4.2. P d denote scattered disturbance features due to footpaths, the black 
arrows indicate high correlation aligned along the disturbance tracks d2 and d5, as showed in Fig. 4, e 1 indicate areas affected by potential erosion processes. 


visible in the 2005 and 2013 images. This is due to a number of factors 
such as the different sensors (see Table 1), time of acquisition, illumi¬ 
nation conditions, noise content and disturbance of tracks. Moreover, 
the 2011 and 2013 images are more noisy respect to the 2002 and 2005 
scenes. To reduce the noise the filtering will inevitably smooth edges of 
the geoglyphs, making them less visible, whereas the use of PCA 
minimizes the noise effect without reducing the spatial resolution. We 
have also to consider that the diverse wind conditions impacts on the 
visibility of the geoglyph edges due to the movement of desert sands. In 
fact, comparing the 2002 and 2005 scenes, acquired by the same sensor 
(QuickBird) and during the same time of the day (about 15:17 and 
15.29, respectively) but, in diverse periods of the year, September and 
March, respectively, the visibility of the geoglyphs is different. In par¬ 
ticular, it is better in March than in September (see Fig. 4a and b, re¬ 
spectively), due to the Paracas dust storms which usually affects the 
desert of Nasca mainly in August and in September (Briceno-Zuluaga 
et al., 2017). All these drawbacks are reduced exploiting the average 
effect of PCI, thus enhancing the geoglyph edges. 

A further emphasis of the geoglyph features is obtained by using 
some spatial texture analyses, among which the dissimilarity provided 
the best results. Fig. 8 shows satellite panchromatic scene of 2005 (8a), 
the first component of PCA (8b), the dissimilarity computed for the 
2005 (8c) and the RGB composition of dissimilarity as obtained for the 
single scenes acquired in 2002, 2005, and 2011 (8d). Moreover two 
zoomed details are showed in Fig. 8e-n. They focus the intersection of 
several geoglyphs. The visual comparison of the 4 scenes and zoomed 
details evidence the added value of dissimilarity respect to both the 
2005 image and PCI in revealing some details which helped us to un¬ 
derstand the relative temporal sequence of geoglyphs. 

Both PCI and dissimilarities computed for the single scenes have 
been classified and followingly segmented in order to automatically 
extract geoglyph features. However, the results obtained were not sa¬ 
tisfactory. In particular, the geoglyphs with a higher discriminability 
(that is with a better state of conservations) have been extracted. 
Whereas other subtle geoglyph features were not extracted. Finally, a 
huge number of false alarms make the classification useless. 

On the basis of the panchromatic satellite imagery (acquired on 
2002, 2005, 2011 and 2013) and derived products (PCI and 
Dissimilarity), a comprehensive map of the geoglyphs has been ob¬ 
tained as displayed in Fig. 9a. Fig. 9b is a subset bordered by black 
rectangle which in Fig. 9a puts in evidence several lacunas of the 
geoglyphs due to disturbance tracks caused by vehicle and motorcycles, 


erosion processes and footpaths. These lacunas have been digitally in¬ 
tegrated allowing us to obtain the original configuration of the geo¬ 
glyphs. Such ‘digital restoration’ (i.e. a virtual reconstruction of geo¬ 
glyphs) could be useful in the future for a real restoration of the 
geoglyphs, object of debate especially for the most famous biomorphic 
geoglyphs. 

The availability of different products, in particular PCI and dis¬ 
similarity, greatly facilitated the reconnaissance of the geoglyphs which 
is not an easy task, even for who had the opportunity to observe them 
on ground during field survey because they can be confused with other 
anthropic traces. In order to reduce potential errors in the re¬ 
connaissance of geoglyphs the investigation strategy also included a 
validation phase performed by close range surveys conducted using 
aircrafts and drones. These enabled a deeper analysis of geoglyphs 
adding details from the higher resolution DEMs and orthophotos. They 
were obtained by SfM-based processing of aerial and drone images (see 
the aerial images by aircraft and drone in Fig. 2d, e, 2f). The area 
surveyed by aircrafts (with spatial resolution of DEM and orthophoto 
equal to 31.7 cm, and 15.9 cm, respectively) was the centre of the 
ground drawing scene, including in particular the meandering motif 
(see Fig. 2ab) selected to cover this particular significant area. Whereas, 
the subset taken from drone (with spatial resolution of DEM and or¬ 
thophoto equal to 5 cm, and 2.5 cm, respectively) was centred on an 
area characterized by the intersection of diverse geoglyphs (see 
Fig. 2ab). The DEMs have been enhanced using SVF and SLRM to em¬ 
phasize edges and microtopography. In particular, SLRM well enhanced 
the borders of the geoglyphs because it facilitates the perception and 
discrimination of sunken from raised features and emphasizes height- 
differences between individual features and the surrounding area. SVF 
mainly delineated concave microrelief preserving the perception of the 
general microtopography. 

In particular, the high detail of drone based orthophotos, DEMs and 
derived models were particularly useful for the spatial characterization 
of the intersection of geoglyphs. Fig. 10 shows a 3d detail of DEM, 
Orthophoto, LRM and SVF centred on the overlapping of three different 
geoglyphs: one trapezoid, indicated with red arrows, and two rectan¬ 
gular figures indicated with white and light blue arrows. The four 
images in Fig. 11a, b, 11c, lid helped us in the identification of the 
temporal sequence (before and after) of the three geoglyphs of the 
subset. In particular, if we focus on the geoglyphs spatial continuity we 
can infer their relative temporal sequence. On the basis of this criterion, 
it is reasonable to assume that the trapezoid, indicated with red arrows. 


12 



























N. Masini and R. Lasaponara 


Remote Sensing of Environment 236 (2020) 111 447 



0 300m 



0 _ 100m 


Fig. 8. (a) Satellite panchromatic scene of 2005, (b) the first component of PCA, (cj dissimilarity computed for the 2005, and Id) the RGB composition of dissimilarity 
as obtained for the single scenes acquired in 2002, 2005, and 2011. The two subsets (from e to n) provide the details for the two subareas related to linear and 
curvilinear geoglyph features in (e) to (h), and (i) to (n) respectively as in the previous figures (a-d). 


has been done before the rectangular figure indicated with light blue 
arrows. Similarly, we can assume that the rectangular element in¬ 
dicated with white arrows has been done after the trapezoid. If we focus 
on the cross between the rectangular figures we can argue that the 
rectangle indicated with light blue arrows was realized after the rec¬ 
tangle indicated with light blue arrows. Therefore, the temporal se¬ 
quence in the execution of the three geoglyphs was the following: Phase 
I trapezoid, phase II rectangle indicated with white arrows, phase III 
rectangle indicated with light blue arrows. 

This approach extended to the entire scene of geoglyphs allowed us 
to have a comprehensive relative chronology of the area under 


investigation (see Fig. 11) which put in evidence four different phases 
of execution of the geoglyphs around the meandering figure. 

As a whole, Fig. 12 shows the methodological approaches we de¬ 
vised for enhancing geoglyphs and for the detection and automatic 
extraction of disturbance features. Moreover, Fig. 12 also provides an 
overview of the outputs from the diverse analyses herein conducted and 
points out the most significant results. 

5. Conclusions and future perspectives 

The availability and open access to archive and novel data from 


13 











N. Masini and R. Lasaponara 


Remote Sensing of Environment 236 (2020) 111 447 



Fig. 9. (a) Map of geoglyphs, (b) ‘Virtual reconstruction’ of geoglyphs showed in Fig. 9a bordered by black rectangular box. 



Fig. 10. 3D visualization of UAV based DEM (a), 
orthophoto (b), and derived texturized models 
obtained by computing SLRM (c) and SVF (d). 
The coloured arrows identify the diverse phase 
of drawing of the geoglyphs. The subset is lo¬ 
cated in Fig. 2d bordered by green box. (For 
interpretation of the references to colour in this 
figure legend, the reader is referred to the Web 
version of this article.) 
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Phase I 



Phase III 



Phase II 



Phase IV 




Fig. 11 . Relative chronology of execution phases of geoglyphs of subset in Fig. 9b. For other geoglyphs it has not been possible to understand diachronic relation with 
other figures. 


spaceborne and airborne platforms are today strongly increasing, thus 
favoring the use, reuse and integration of EO data in the cultural 
heritage (CH) domain, including archaeological heritage from the dis¬ 
covery to the monitoring, from the enhancement to the preservation. 
This requires efforts addressing to manage and process big data, make 
available reliable automatic data processing approaches and integrated 
methodologies of data analysis and interpretation also via case study 
demonstration to bridge EO technology and archaeological commu¬ 
nities. 

This paper deals with a multiscale, multisensor and multitemporal 
RS based approach for the analysis, interpretation and monitoring of 
Nasca geoglyphs which are one of the most impressive and vulnerable 
examples of cultural heritage throughout the entire world. They are 
exposed to damage much more than any other cultural heritage, due to 
their intrinsic fragility but also because they are characterized by very 
subtle features and strongly threatened by several anthropogenic fac¬ 
tors. 

The area investigated in this paper is located in Pampa de Atarco, 
selected because it is representative of the most fragile geoglyphs, much 
more than those already investigated in other sites, such as Pampa de 
Nasca (Hesse, 2015; Comer, 2017) and Palpa (Lambers, 2006). This is 
due to the specific engraving techniques based on the use of sand grain, 
pebbles and gravel of small dimension which make the geoglyphs very 
difficult to be detected not only on ground but also by satellite and close 
range view, especially in not optimal illumination condition. This has 
suggested an investigation strategy based on the integration of multi¬ 
scale remote sensing passive and active data, including close range 
imagery, to enhance the visibility of geoglyphs (most of them are very 
subtle) and to capture the changes due to human induced disturbance 
affecting them. To cope with these issues, the analyses herein per¬ 
formed were based on statistical indicators computed in both temporal 
and spatial domain. Multitemporal analysis of VHR satellite imagery 


using PCA and Skewness allowed us to enhance the visibility of dis¬ 
turbance features and to automatically extract them using unsupervised 
classifications. In particular, respect to the satellite multitemporal data 
set which reveals the presence of five off-road vehicle tracks, the en¬ 
hancement by PCA and Skewness enabled to add other five thinner 
tracks referable to footpaths or motorcycle tracks. The best results in 
terms of enhancement and automatic extraction capability of dis¬ 
turbance features have been obtained by Skewness. 

The reuse of UAV L SAR-based correlation map, available free of 
charge from NASA, provided useful information on the state of dis¬ 
turbance from 2013 to 2015, widening the observation time window of 
the VHR satellite data set from 2002 to 2013. 

In the spatial domain, the textural analyses mainly based on dis¬ 
similarity along with PCI enabled us to enhance edges of geoglyphs 
features facilitating their mapping and interpretation. In particular, in 
areas characterized by the superposition of geoglyphs, thanks also to 
the contribution of close range UAV photogrammetry, it was possible to 
reconstruct the relative chronological sequence thus providing an im¬ 
portant contribution to archaeological studies. In the future, it is de¬ 
sirable to improve the relative chronological method with approaches 
and tools for absolute dating. This can be achieved by the analysis of 
pottery and ceramics close and along the geoglyphs, previously iden¬ 
tified using Remote sensing data. 

As a whole, results from our investigations pointed out that the 
integration of satellite-based imaging and information with manned 
and unmanned data have significant potential to contribute into deci¬ 
sion making processes for improving knowledge and supporting sys¬ 
tematic monitoring and sustainable management strategies as manda¬ 
tory steps to preserve unique or un-renewable resources as cultural 
heritage and landscape. 
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Graphical summary of methodological approach and outputs 
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Fig. 12. Graphic summary of the two multiscale, multisensor and multitemporal approaches: one devised for enhancing geoglyphs, the second for detecting and 
automatically extracting disturbance features. 
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APPENDIX 

Principal Component Analysis 

The PCA is based on the covariance (or variance) matrix (S) computed among all input time series, the eigenvalues, eigenvectors and loadings 
using formulas la-4 


covkik 2 = 


n x m . 


X X (SBi,j, kl P k2 ) x (SBjj ik2 p k2 ) 


i=i j=i 


[la] 


where k 1; k 2 are two input spectral channels, SBij is the spectral value of the given channel in row i and column), n number of row, m number of 
columns, p is the mean of all pixel SBy values in k 2 and k 2 

The most common application of PCA is based on the use of covariance matrix, even if both variance and covariance measure the “spread” of a 
dataset around the mean. A series of new image layers (eigenvectors generally called eigenchannels or components) are obtained from the eigen¬ 
values of the covariance matrix (using formula 2 and 3). The percentage of total dataset variance, explained by each component, is obtained by 

formula 2: 


% i = 100 x Sj/ Yj h 


[ 2 ] 


where % ij are eigenvalues of S. 

Finally, a series of new image layers (called eigenchannels or components) are computed (formula 3) by multiplying, for each pixel, the ei¬ 
genvector of S for the original value of a given pixel in the input bands 


p i = X Pk x U W 


[3] 


where Pi indicates a spectral channel in component i, u ki eigenvector element for component i in input band k, P k spectral value for channel k, 
number of input band. 

A loading, or correlation R, of each component i with each input date k can be calculated by using formula 4. 


Rkj = u k ,i x (A ;)2 x (var k ) 1/2 

where var k is the variance of input date k (obtained by reading the kth diagonal of the covariance matrix. 

Skewness measurement 

The most common skewness measurement is Fisher-Pearson (or Pearson median skewness) expressed by formula 5 
Skewness = 


[4] 


w U 


Standard deviation 


[5] 


where Xj is the value of the pixel for band j, N is the number of the images, in our case satellite scenes acquired in a different time, xi is the mean 

Standard deviation = (variance) 1 ^ 

N 


Variance = 


= —-— X ( x i — x i ) 2 

N-l 1 


[ 6 ] 


Texture Indicators 

The Variance that measures local variance of the processing window 
Variance = ^ (i - u) 2 (pi, j ) 


• The Correlation that measures the linear dependency of a gray level on those of neighbouring image cells. 
Z;ZjP('J)- m x my 


Correlation = 


s x sy 


• Contrast that measures gray-level contrast by using weighting factors equal to the square of the gray level difference. 


[7] 


[ 8 ] 
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Ng-l Ng Ng 

Contrast = Z (rc) 2 ZZ P(ij ) 

„=o i=i j=i [9] 

• Dissimilarity that uses weighting factors equal to the absolute value of the gray level difference, so 

Ng-l Ng Ng 

Dissimilarity = Z (n) EE P(U) 2 

n =0 i =1 j=l [10] 

• Homogeneity measures the smoothness of the image area. 

Homogeneity = ± + J _ J)2 p(i, j) 

• Entropy that measures the randomness or disorder of the image area. 

Entropy = Z ^p(i,7)log(py)) 

i j=i [12] 

Where, in equations 7 to 12, i, j are coordinates of the co-occurrence matrix space; p (i, j) is the element of the i and j coordinates; and Ng is 
dimension of the co-occurrence matrix, which has gray value range of the original image. Finally, in equation (8) m and s denote average and 
standard deviation, respectively. 
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